Gas exchange in Delmarva bays

library(CHNOSZ)
data(thermo)
library(xts)
library(dplyr)
library(dygraphs)
library(purrr)
library(ggplot2)
library(lubridate)
library(readr)
library(data.table)
library(forcats)

Calculate a time series of \(k\) in \(\frac{cm}{hr}\) using inputs: water temperature, water depth, wind speed. \(k\) is calculated using equation 2 in Poindexter et al. 2016, which uses a different formulation for wind-driven vs. thermal convection driven gas exchange based on thresholds for heat flux.

Background

Gas flux equation

\(Flux = k(C_{eq} - C_w)\)

\(k\) is typically reported as \(k_{600}\) to facilitate comparison.

\(k_{600} = k(\frac{600}{Sc})^{-n}\)

\(n\) is a factor characterizing the kinematic behavior of the water surface.

Previous studies

The most relevant studies on calculating gas exchange are:

• Read et al. 2012
• Holgerson et al 2017
• MacIntyre et al. 2010 - model \(k\) with two functions of \(U_{10}\), one for periods of surface heat loss and one for periods of heat gain
• Soloviev et al. 2007 - combined effects of surface heat loss and wind via surface dissipation of turbulent kinetic energy (from combined effects of buoyancy, wind shear, wave breaking)
• Poindexter and Variano 2013 - develop an analytical relationship between surface heat flux \(q\) and \(k\) (based on laboratory measurements)
• Poindexter et al. 2016

Poindexter and Variano (2013) developed new models for wetlands with emergent vegetation, based on the idea that emergent vegetation would attenuate wind speed above the water surface, modify fluid shear at the water surface, and influence stirring beneath the water surface. They focus on thermal convection, parameterized from surface heat loss, and wind shear, parameterized using the mean wind speed in the constant-velocity or shear-free region of the vegetation canopy.

Forested wetlands like Delmarva bays can be inundated areas without emergent vegetation (e.g. forested ponds), where tree canopy attenuates wind speed but effects of plant stems are not as relevant.

Using Poindexter et al 2016 method for wetland perspective, where emergent vegetation and surfactants are important. New relationships for \(k\) designed specifically for use in the presence of an emergent vegetation canopy. Based on lab measurements with a model wetland, which suggested a quadratic relationship between \(k\) and the mean in-canopy wind speed (Poindexter and Variano, 2013).

\(k\) is wind-driven when \(q \geq 0\) or \(Ra < 8 \times 10^6\), and thermal convection driven when \(q < 0\) and \(Ra > 8 \times 10^6\). \(Ra\) is the Rayleigh number, a function of heat flux. When \(Ra\) is below a critical value for a fluid, then heat transfer is primarily in the form of conduction; when it exceeds that value then heat transfer is primarily in the form of convection. *Except if water temperature is below 4 deg C, then \(\beta\) changes sign and \(k\) is thermally driven by \(q > 0\).

• First check if \(q\geq 0\), if yes, use \(k_{wind}\) unless water temp < 4.
• If \(q < 0\), use \(k_{thermal}\) unless \(Ra < 8\times 10^6\).

Surface Heat Flux

Surface heat flux \((q)\) can be used to estimate thermal convection. Negative heat fluxes indicate a cooling water column (water losing heat to the air) and positive heat flux means the water column is warming. Typical values of \(q\) in temperate wetlands are -200 to 300 W m^-2. Convective mixing occurs when the water is losing heat to the air \((q < 0)\).

\(q = \frac{dT_b}{dt}\rho c_pH\)

• \(c_p\) is the isobaric heat capacity (isobaric = constant pressure), in \(\frac{J}{kgK}\)
• \(\rho\) is water density, in \(\frac{kg}{m^3}\)
• \(H\) is the depth of the water column, in meters
• \(T_b\) is the bulk water temperature

Calculate \(q\), using functions from the CHNOSZ package to calculate water density and isobaric heat capacity

calc_q_heatflux <- function(temp_C, waterdepth_m, dT, dt_sec = 3600){
  # isobaric heat capacity cal K^-1 mol^-1
  cp = as.numeric(
    CHNOSZ::water("Cp", T = 273.15 + temp_C))
  # water density kg m^-3
  rho = as.numeric(
    CHNOSZ::water("rho", T = 273.15 + temp_C))
  # convert cp to J kg^-1 K^-1
  cp = cp*(4.184/0.018)
  
  q_wm2 <- (dT/dt_sec) * rho * cp * waterdepth_m
  return(q_wm2)
}

Rayleigh Number calculation

Used in combination with \(q\) to determine which formula to use for calculation of \(k\).

\(Ra = 4.7(\frac{-qg\beta L^4}{\alpha^2vc_p\rho})^{\frac{3}{4}}\)

• \(g\) - acceleration due to gravity
• \(q\) - water column heat flux
• \(L\) - characteristic length scale, wetland depth

Poindexter and Variano (2013) find that using the wetland depth as the length scale \(L\), the Rayleigh number threshold is met for wetland water columns deeper than 10cm for typical heat loss rates.

• \(\beta\) - thermal expansion coefficient
• \(\alpha\) - thermal diffusivity
• \(v\) - kinematic viscosity
• \(c_p\) - heat capacity at constant pressure
• \(\rho\) - density

calc_rayleigh <- function(q, temp_C = 25, lscale = 1){
  g = -9.8 # m s^-2
  beta = as.numeric(
    CHNOSZ::water("alpha", T = 273.15 + temp_C)) # K^-1
  alpha = as.numeric(
    CHNOSZ::water("tdiff", T = 273.15 + temp_C)) #cm^2 s^-1
  kv = as.numeric(
    CHNOSZ::water("visck", T = 273.15 + temp_C)) #cm^2 s^-1
  cp = as.numeric(
    CHNOSZ::water("Cp", T = 273.15 + temp_C))
  rho = as.numeric(
    CHNOSZ::water("rho", T = 273.15 + temp_C)) # kg m^-3
  # convert cm to m
  alpha = alpha/10000
  kv = kv/10000
  
  ra <- 4.7*((-q*g*beta*lscale^4)/
       (alpha^2*kv*cp*rho))^(3/4)
  return(ra)
}

calc_rayleigh(q = .1, temp_C = 30, lscale = .10)

[1] 4532289

calc_rayleigh(q = .10, temp_C = 30, lscale = .10) > 8e6

[1] FALSE

Thermal convection K

From Poindexter et al. 2016 (GRL), a semi-emprical function derived from the relationship between heat flux and thermal convection, which is similar to the surface renewal models used for periods of relative calm in lakes and oceans by Read et al. 2012 and Soloviev et al. 2007.

\(k = 0.14^{\frac{3}{4}}(\frac{-qg\beta\alpha^2}{vc_p\rho})^{\frac{1}{4}}(\frac{Sc}{Pr})^{-n}\)

• \(q\) - water column heat flux, calculated from rate of change in wetland water column temperature with time (described in supporting info)
• \(g\) - acceleration due to gravity
• \(Sc\) - Schmidt number

Functions of water temperature and atmospheric pressure:

• \(\beta\) - thermal expansion coefficient
• \(\alpha\) - thermal diffusivity
• \(v\) - kinematic viscosity
• \(c_p\) - heat capacity at constant pressure
• \(\rho\) - density
• \(Pr\) - Prandtl number

Calculate k for thermal convection

calc_k_poindexter2b <- function(q, temp_C, n_surface = (2/3)){
  # if(q < 0){stop("Negative heat flux. 
  #                Please use wind model to calculate K.")}
  g = 9.8 # m s^-2
  beta = as.numeric(
    CHNOSZ::water("alpha", T = 273.15 + temp_C)) # K^-1
  alpha = as.numeric(
    CHNOSZ::water("tdiff", T = 273.15 + temp_C)) #cm^2 s^-1
  kv = as.numeric(
    CHNOSZ::water("visck", T = 273.15 + temp_C)) #cm^2 s^-1
  cp = as.numeric(
    CHNOSZ::water("Cp", T = 273.15 + temp_C))
  rho = as.numeric(
    CHNOSZ::water("rho", T = 273.15 + temp_C)) # kg m^-3
  Pr = as.numeric(
    CHNOSZ::water("Prndtl", T = 273.15 + temp_C))
  Sc = 600
  
  # convert cp to J kg^-1 K^-1
  cp = cp*(4.184/0.018)
  # convert cm to m
  alpha = alpha/10000
  kv = kv/10000
  
  k = (0.14^(3/4)) *
    ((-q * g * beta * alpha^2)/(kv * cp * rho))^(1/4) * 
    (Sc/Pr)^(-n_surface)
  
  k_cmhr = k*100*3600
  return(k_cmhr)
}

calc_k_poindexter2b(q = 50, temp_C = 25) # cm/hr

[1] NaN

calc_k_poindexter2b(q = -50, temp_C = 25) # cm/hr

[1] 0.636363

Wind-driven k

\(k = 3\Gamma{\langle U_{canopy}}\rangle^2 (\frac{Sc}{600})^{-n}\)

• \(\Gamma\) enhancement factor to account for fluctuations that increase the gas transfer velocity relative to value predicted from mean wind speed alone. Poindexter et al. use \(\Gamma = 0.32\) based on in-canopy wind speed measurements at their marsh study site. (Poindexter and Variano use 1.3 though?)
• \(U_{canopy}\) is mean in-canopy wind speed, which can be computed from wind shear stress and wind speed measured above the canopy (but is being monitored at sites QB and ND)

Calculate k for wind driven turbulence.

calc_k_poindexter2a <- function(ws, gamma = (1/3), Sc = 600, n_surface = (2/3)){
  k = 3*gamma*(ws)^2 * (Sc/600)^(-n_surface)
  return(k)
}

In-canopy wind speed is assumed to be in the zone where wind speed profile is vertical?

Data

Load data to use the equations above for calculating \(k\). Aggregate high resolution data to hourly means.

Temperature

Water temperature measured at 3 locations in a vertical profile, in 2 wetlands. Need to compare with water level to determine when sensors were in the water and when they were exposed (also evident by looking at daily fluctuations).

Look at typical daily time series pattern of temperature in the water column, to determine how to calculate \(\frac{dT}{dt}\). Poindexter and Variano use linear regression of temperature readings during their experiments.

temp_df <- read_csv("data/watertemp_qbMID.csv")

Parsed with column specification:
cols(
  datetime = col_datetime(format = ""),
  temp_C = col_double()
)

Wind

In-canopy wind speed measured with sonic anemometers approximately 1m over the water surface (varies between xx and xx due to changes in water level).

wind_df <- read_csv("data/windspeed_qb.csv")

Parsed with column specification:
cols(
  datetime = col_datetime(format = ""),
  ws_ms = col_double()
)

wind_df %>% 
  mutate(date = lubridate::date(datetime)) %>%
  mutate(time_hr = lubridate::hour(datetime)) %>% 
  # filter(date >= as.Date('2018-06-01') & date <= as.Date('2018-07-01')) %>%
  ggplot(aes(x = time_hr, y = ws_ms, group = date)) +
  geom_line(aes(col = date)) +
  # facet_wrap(~ lubridate::month(date)) +
  # facet_wrap(~ lubridate::yday(date)) +
  # facet_wrap(~ lubridate::wday(date, TRUE)) +
  xlim(c(0, 24)) +
  theme_bw()

wind-driven k is just a function of wind speed

wind_df <- wind_df %>%
  mutate(k_wind = map_dbl(ws_ms, ~calc_k_poindexter2a(ws = .x)))

wind_df %>%
  mutate(date = lubridate::date(datetime)) %>%
  mutate(time_hr = lubridate::hour(datetime)) %>% 
  # filter(date >= as.Date('2018-06-01') & date <= as.Date('2018-07-01')) %>%
  ggplot(aes(x = time_hr, y = k_wind, group = date)) +
  geom_line(aes(col = date)) +
  # facet_wrap(~ lubridate::month(date)) +
  # facet_wrap(~ lubridate::yday(date)) +
  # facet_wrap(~ lubridate::wday(date, TRUE)) +
  xlim(c(0, 24)) +
  theme_bw()

wind_df %>%
  mutate(date = lubridate::date(datetime)) %>%
  mutate(time_hr = lubridate::hour(datetime)) %>% 
  mutate(wind_thresh = ws_ms > 0.7) %>%
  # filter(date >= as.Date('2018-06-01') & date <= as.Date('2018-07-01')) %>%
  ggplot(aes(x = datetime, y = k_wind, group = date)) +
  geom_point(aes(col = wind_thresh), size = 0.5) +
  # facet_wrap(~ lubridate::month(date)) +
  # facet_wrap(~ lubridate::yday(date)) +
  # facet_wrap(~ lubridate::wday(date, TRUE)) +
  # xlim(c(0, 24)) +
  theme_bw()

Water depth

wl_df <- read_csv("data/waterlevel_qb.csv")

Parsed with column specification:
cols(
  date = col_date(format = ""),
  wl_m = col_double(),
  datetime = col_datetime(format = "")
)

Combine data

Combine 3 data sources into one data frame (wl_df, wind_qb, temp_df)

head(wl_df) # water level is daily

# A tibble: 6 x 3
  date        wl_m datetime           
  <date>     <dbl> <dttm>             
1 2017-10-01 0.629 2017-10-01 04:00:00
2 2017-10-02 0.619 2017-10-02 04:00:00
3 2017-10-03 0.610 2017-10-03 04:00:00
4 2017-10-04 0.603 2017-10-04 04:00:00
5 2017-10-05 0.593 2017-10-05 04:00:00
6 2017-10-06 0.583 2017-10-06 04:00:00

head(wind_df) # hourly wind speed starts 6-27, includes k

# A tibble: 6 x 3
  datetime            ws_ms k_wind
  <dttm>              <dbl>  <dbl>
1 2018-06-27 13:55:00 0.580  0.336
2 2018-06-27 14:55:00 0.616  0.379
3 2018-06-27 15:55:00 0.549  0.302
4 2018-06-27 16:55:00 0.608  0.370
5 2018-06-27 17:55:00 0.747  0.558
6 2018-06-27 18:55:00 0.685  0.469

head(temp_df) # hourly temp, starts 7-21

# A tibble: 6 x 2
  datetime            temp_C
  <dttm>               <dbl>
1 2018-07-21 17:58:44   26.9
2 2018-07-21 18:58:44   25.7
3 2018-07-21 19:58:44   25.2
4 2018-07-21 20:58:44   24.7
5 2018-07-21 21:58:44   24.2
6 2018-07-21 22:58:44   23.8

tail(temp_df) # ends 9-30

# A tibble: 6 x 2
  datetime            temp_C
  <dttm>               <dbl>
1 2018-09-30 09:58:44   18.0
2 2018-09-30 10:58:44   17.8
3 2018-09-30 11:58:44   17.7
4 2018-09-30 12:58:44   17.5
5 2018-09-30 13:58:44   17.5
6 2018-09-30 14:49:59   NA  

tail(wind_df) # ends 9-30

# A tibble: 6 x 3
  datetime            ws_ms k_wind
  <dttm>              <dbl>  <dbl>
1 2018-09-30 11:55:00 0.234 0.0546
2 2018-09-30 12:55:00 0.240 0.0574
3 2018-09-30 13:55:00 0.25  0.0625
4 2018-09-30 14:55:00 0.302 0.0910
5 2018-09-30 15:55:00 0.273 0.0747
6 2018-09-30 16:20:00 0.252 0.0635

tail(wl_df) # ends 9-04

# A tibble: 6 x 3
  date        wl_m datetime           
  <date>     <dbl> <dttm>             
1 2018-08-30 0.380 2018-08-30 04:00:00
2 2018-08-31 0.369 2018-08-31 04:00:00
3 2018-09-01 0.361 2018-09-01 04:00:00
4 2018-09-02 0.350 2018-09-02 04:00:00
5 2018-09-03 0.340 2018-09-03 04:00:00
6 2018-09-04 0.328 2018-09-04 04:00:00

hourly time series of wind, water level, temp between 7-21 and 9-04

Rolling merge with data.table package

temp_dt <- temp_df %>% 
  data.table() %>%
  setkey(datetime)
wind_dt <- wind_df %>% 
  data.table() %>%
  setkey(datetime)
kdata <- temp_dt[wind_dt, roll = "nearest"]
kdata <- kdata %>%
  mutate(date = lubridate::date(datetime)) %>% 
  left_join(dplyr::select(wl_df, -datetime))

Joining, by = "date"

write_csv(kdata, "data/kdata.csv")

Apply Calculations

Make a function that will use that as input to calculate a new column with \(k\) values.

1. Calculate dT from temperature column
2. Add temp_thresh (if temp_C < 4)
3. Calculate q (heat flux) and q_neg (if \(q < 0\))
4. Calculate Ra and Ra_thresh (if \(Ra > 8\times 10^6\))
5. Determine k_category (wind or thermal), if conditions for \(k_{thermal}\) are met: q_neg = TRUE, Ra_thresh = TRUE, or temp_thresh = TRUE.
6. For rows where k_category == "wind", calculate k_cmhr using calc_k_poindexter2a().
7. For rows where k_category == "thermal", calculate k_cmhr using calc_k_poindexter2b().

Calculate time series of \(q\) surface heat flux (\(Wm^{-2}\)) as a function of temperature time series and water column depth

• \(\frac{dT_b}{dt}\) in \(\frac{K}{sec}\) will need to convert from \(\frac{C}{hr}\)

kdf <- kdata %>% 
  filter(!is.na(wl_m)) %>%
  # mutate(q_neg = q_wm2 <0) %>%
  mutate(temp_thresh = temp_C < 4)

head(kdf)

             datetime   temp_C     ws_ms    k_wind       date      wl_m temp_thresh
1 2018-06-27 13:55:00 26.91349 0.5800000 0.3364000 2018-06-27 0.5550974       FALSE
2 2018-06-27 14:55:00 26.91349 0.6158333 0.3792507 2018-06-27 0.5550974       FALSE
3 2018-06-27 15:55:00 26.91349 0.5491667 0.3015840 2018-06-27 0.5550974       FALSE
4 2018-06-27 16:55:00 26.91349 0.6083333 0.3700694 2018-06-27 0.5550974       FALSE
5 2018-06-27 17:55:00 26.91349 0.7466667 0.5575111 2018-06-27 0.5550974       FALSE
6 2018-06-27 18:55:00 26.91349 0.6850000 0.4692250 2018-06-27 0.5550974       FALSE

Crop to just July 24 - August 24

kdf <- kdf %>% 
  dplyr::filter(date >= as.Date('2018-07-23') & date <= as.Date('2018-08-24'))

kdf %>% 
  dplyr::select(-date, -temp_thresh, -k_wind) %>%
  tidyr::gather(variable, value, temp_C:wl_m) %>%
  ggplot(aes(x = datetime, y = value)) +
  geom_line() +
  facet_wrap(vars(variable), scales = "free_y", ncol = 1) +
  theme_bw()

Heat flux

kdf$dtemp <- c(NA, diff(kdf$temp_C)) # revisit for smoothing
kdf %>%
  # slice(-1) %>% 
  mutate(dt_pos = dtemp>0) %>%
  ggplot(aes(x = datetime, y = temp_C)) +
  geom_line(col = "gray") +
  geom_point(aes(col = dt_pos)) + theme_bw()

kdf$q_wm2 <- pmap_dbl(list(temp_C = kdf$temp_C,
              waterdepth_m = kdf$wl_m,
              dT = kdf$dtemp),
              calc_q_heatflux)

kdf %>%
  mutate(dt_pos = dtemp>0) %>%
  ggplot(aes(x = datetime, y = q_wm2)) +
  geom_line(col = "gray") +
  geom_point(aes(col = dt_pos)) + theme_bw()

# kdf$Ra <- kdf %>% pmap_dbl(calc_rayleigh())
kdf$Ra <- pmap_dbl(list(q = kdf$q_wm2, 
                        temp_C = kdf$temp_C, 
                        lscale = kdf$wl_m), 
                   calc_rayleigh)

kdf$Ra_thresh <- kdf$Ra > 8e6

kdf %>%
  ggplot(aes(x = datetime, y = Ra)) +
  geom_line(col = "gray") +
  geom_point(aes(col = Ra_thresh)) + theme_bw()

# kdf <- 
  kdf %>%
  mutate(q_pos = q >=0) %>% View()

Error in mutate_impl(.data, dots) : 
  Evaluation error: comparison (5) is possible only for atomic and list types.

purrr::pmap_dbl(list(q = kdf_thermal$q_wm2, temp_C = kdf_thermal$temp_C),
                   calc_k_poindexter2b) %>% head()

[1] 0.2875823 0.3616141 0.3703332 0.3968111 0.3660502 0.3744688

kdf_join %>%
  mutate(k_mix = dplyr::case_when(
    k_category == "wind" ~ k_wind,
    k_category == "thermal" ~ k_thermal
  )) %>% head(10)

              datetime   temp_C     ws_ms    k_wind       date      wl_m temp_thresh        dtemp     q_wm2         Ra
1  2018-07-23 00:55:00 22.55509 0.3558333 0.1266174 2018-07-23 0.4849155       FALSE           NA        NA         NA
2  2018-07-23 01:55:00 22.55139 0.3812500 0.1453516 2018-07-23 0.4849155       FALSE -0.003703704 -2.083821        NaN
3  2018-07-23 02:55:00 22.54213 0.4425000 0.1958063 2018-07-23 0.4849155       FALSE -0.009259259 -5.209564        NaN
4  2018-07-23 03:55:00 22.54815 0.4829167 0.2332085 2018-07-23 0.4849155       FALSE  0.006018519  3.386212 5449991329
5  2018-07-23 04:55:00 22.53796 0.6050000 0.3660250 2018-07-23 0.4849155       FALSE -0.010185185 -5.730526        NaN
6  2018-07-23 05:55:00 22.52454 0.5762500 0.3320641 2018-07-23 0.4849155       FALSE -0.013425926 -7.553898        NaN
7  2018-07-23 06:55:00 22.53472 0.5154167 0.2656543 2018-07-23 0.4849155       FALSE  0.010185185  5.730530 8081523609
8  2018-07-23 07:55:00 22.54167 0.5283333 0.2791361 2018-07-23 0.4849155       FALSE  0.006944444  3.907173 6065690561
9  2018-07-23 08:55:00 22.53194 0.5337500 0.2848891 2018-07-23 0.4849155       FALSE -0.009722222 -5.470055        NaN
10 2018-07-23 09:55:00 22.52130 0.5654167 0.3196960 2018-07-23 0.4849155       FALSE -0.010648148 -5.991027        NaN
   Ra_thresh k_category q_pos k_thermal     k_mix
1         NA       <NA>    NA        NA        NA
2         NA    thermal FALSE 0.2875823 0.2875823
3         NA    thermal FALSE 0.3616141 0.3616141
4       TRUE       wind  TRUE        NA 0.2332085
5         NA    thermal FALSE 0.3703332 0.3703332
6         NA    thermal FALSE 0.3968111 0.3968111
7       TRUE       wind  TRUE        NA 0.2656543
8       TRUE       wind  TRUE        NA 0.2791361
9         NA    thermal FALSE 0.3660502 0.3660502
10        NA    thermal FALSE 0.3744688 0.3744688

Notes

Ideas

Compare methane flux and k measurements with and without changing water levels - what is the effect of including changing water level? what is the relative importance of a) changing areal extent of water surface to scale flux by? b) change in k value and k category on the magnitude of the flux? c) dilution effect of the methane concentration in the water?

Surface heat flux approach in Read et al. 2012 assumes that sediment heating, advected components, and changes in water levels are negligible. These assumptions do not hold in delmarva bay wetland systems!

Delmarva bay wetland sites are sheltered by tree canopy which affects the wind shear driven gas exchange. Read et al. method using approach of scaling the effect of wind sheltering where the size of the lake is small enough, so k is reduced for area affected by sheltering as a function of distance from tree canopy. Unlike that method, there is no water surface area that is unaffected by effects of canopy sheltering. Canopy height measurements from Dubayah et al. CMS product.

In-canopy wind speed is calculated in Poindexter et al. method, as a function of wind speed measured at 10m.

Our measurements of wind speed can be used to evaluate the other methods, by comparing measured in-canopy wind speed to wind speed from a nearby weather station (measured above the canopy). Need to find a local weather station to do that?

Need to understand what the enhancement factor is all about - can I just use \(k = \langle U_{canopy}\rangle^2(\frac{Sc}{600})^{-n}\) - additional factor beyond mean wind speed to characterize effect of wind on \(k\) because the largest scale eddies in the canopy wind field are highly intermittent. Poindexter and Variano 2013 give a relationship for the \(k_{600}\) enhancement factor as a fucntion of \(\frac{\langle(U_{canopy}-\langle U_{canopy}\rangle)^2\rangle^{\frac{1}{2}}}{\langle U_{canopy}\rangle}\). The effect of the gusts is interactions between the wind and plant canopy that generates intermittent, large-scale eddies, which presumably are the gusts that enhance gas exchange in natural wetlands.

Compare wind gust speeeds and wind speeds from sonic anemometer measurements to help think about this.

gusts <- read_csv("data/windgusts_qb.csv")

gusts %>% 
  mutate(date = lubridate::date(datetime)) %>%
  mutate(time_hr = lubridate::hour(datetime)) %>% 
  # filter(date >= as.Date('2018-06-01') & date <= as.Date('2018-07-01')) %>%
  ggplot(aes(x = time_hr, y = gusts_ms, group = date)) +
  geom_line(aes(col = date)) +
  # facet_wrap(~ lubridate::month(date)) +
  # facet_wrap(~ lubridate::yday(date)) +
  # facet_wrap(~ lubridate::wday(date, TRUE)) +
  xlim(c(0, 24)) +
  theme_bw()

Other relevant studies

• Godwin et al. 2013 JGR Evening methane emission events - found that timing of nighttime emission events (as measured by flux towers) coincided with cooling and convective mixing within hollows, which occurred regularly during the growing season. They propose that diurnal thermal stratification in shallow pools traps methane by limiting turbulent transport, and this methane stored during daytime heating is later released during evening cooling due to convective turbulent mixing.
• Happell et al. 1995 used floating chamber measurements in a Florida hardwood swamp and found \(k_{600}\) values averaged \(0.78\pm 0.54 cm h^{-1}\)
• Variano et al. 2009 - \(SF_6\) tracer released in the patterned marshes of the Florida Everglades indicated \(k_{600}\) there ranged from 0.3 to 1.8 cm h^-1
• Sebacher et al. 1983 measured gas exchange in a wetland pond free of emergent vegetation and wind speed 2 cm above the water surface, wind speed around 0.8 m s^-1.